home *** CD-ROM | disk | FTP | other *** search
/ OpenGL Superbible (2nd Edition) / OpenGL SuperBible e2.iso / tools / GLUT-3.7 / PROGS / ADVANCED / VOX.C < prev    next >
Encoding:
C/C++ Source or Header  |  1998-08-12  |  20.5 KB  |  893 lines
  1. /**
  2.  *
  3.  * vox.c : GLUT example for volume rendering
  4.  *
  5.  * Author : Yusuf Attarwala
  6.  *          SGI - Applications
  7.  *
  8.  * Mods by Mark Kilgard.
  9.  *
  10.  * cc vox.c -o vox -lglut -lGLU -lGL -lXmu -lXext -lX11 -lm
  11.  *
  12.  * Voxel Head is a simple volume rendering example using OpenGL 3d textures.
  13.  * This version has limited features (intentional), to keep the code simple.
  14.  *
  15.  * i.e. - there is no control for changing alpha values.
  16.  *      - it just deals with one block of texture (no tiling)
  17.  *      - there are no clipping planes
  18.  *      - it only works with one data set, so it assumes that texture data
  19.  *       is in powers of 2.
  20.  *
  21.  * Some of these features will be implemented in future releases.
  22.  *
  23.  * The technique to create polygonal slices thru voxel space is as
  24.  * follows:
  25.  *
  26.  * Instead of recomputing polygonal slices perpendicular to every
  27.  * viewing vector, it uses the same set of polygonal slices for a
  28.  * range of -45 to 45 degrees. Outside this range, it recomputes
  29.  * another set of slices and so on.
  30.  *
  31.  * These slices are then rendered back to front with 3d texture and
  32.  * blending enabled.
  33.  *
  34.  * It uses GLUT for handling events, windows etc.
  35.  *
  36.  * This program runs with good performance on RealityEngine, VTX,
  37.  * InfiniteReality, or Maximum IMPACT.  In IRIX 6.2, a software
  38.  * implementation of OpenGL's 3D texture mapping extension is
  39.  * available, but performance is very poor.  If you are on a slow
  40.  * machine, you can run "vox -sb" and still get a good feel for how
  41.  * the program performs volume rendering since you can see the 3D
  42.  * slices be rendered back to front with blending.
  43.  * 
  44.  */
  45.  
  46. #include <stdio.h>
  47. #include <stdlib.h>
  48. #ifndef _WIN32
  49. #include <unistd.h>
  50. #else
  51. #define R_OK 04  /* Win32 doesn't define this */
  52. /* Win32 defines these, but with a leading _ */
  53. #define pclose _pclose
  54. #define popen _popen
  55. /* Win32 doesn't have a re-entrand rand() */
  56. #define rand_r rand
  57. #endif
  58. #include <string.h>
  59. #include <math.h>
  60.  
  61. #include <GL/glut.h>
  62.  
  63. #define ABS(a)   (((a) >= 0) ? (a) : -(a))
  64. #define _XY    1
  65. #define _YZ    2
  66. #define _ZX    3
  67. #define _MXY   4
  68. #define _MYZ   5
  69. #define _MZX   6
  70.  
  71. /* pop up menu entries */
  72. #define SPIN_ON   1
  73. #define SPIN_OFF  2
  74. #define MENU_HELP 3
  75. #define MENU_EXIT 4
  76.  
  77. /* global variables */
  78.  
  79. int width, height;      /* window width, height */
  80. float left, right, bottom, top, nnear, ffar;  /* ortho, view volume */
  81. float vol_width, vol_height, vol_depth;  /* volume dimensions */
  82. float bminx, bmaxx, bminy, bmaxy, bminz, bmaxz, bdiag;  /* bounding box */
  83. int n_slices;           /* number of slices */
  84. int tex3dSupported = 1;
  85. float slice_poly[3][4][3];
  86. float slice_tcoord[3][4][3];
  87. float anglex, angley, anglez;
  88.  
  89. unsigned char *voxels;
  90.  
  91. void
  92. readVoxelData(void)
  93. {
  94.   FILE *file;
  95.   int i, j, k, using_pipe;
  96.   unsigned char *vptr;
  97.   unsigned char *vp, *vx;
  98.  
  99.   /* see if the hardware supports 3d texture */
  100. #ifdef GL_EXT_texture3D
  101.   if (!glutExtensionSupported("GL_EXT_texture3D")) {
  102.     printf("\n==================================================================\n");
  103.     printf("This hardware (%s) does not support 3d texture extentions\n",
  104.       (char *) (glGetString(GL_RENDERER)));
  105.     printf("==================================================================\n");
  106.  
  107.     tex3dSupported = 0;
  108.   }
  109. #else
  110.   printf("\n==================================================================\n");
  111.   printf("Not API support for GL_EXT_texture3D extension when compiled.\n");
  112.   printf("==================================================================\n");
  113. #endif
  114.   /* open vox.bin data file */
  115.   if ((file = fopen("vox.bin", "r")) == NULL) {
  116. #ifndef _WIN32
  117.     if (!access("vox.bin.gz", R_OK)) {
  118.       if ((file = popen("gzcat vox.bin.gz", "r")) == NULL) {
  119.         fprintf(stderr, "cannot popen input file vox.bin.gz (missing gzcat?)\n");
  120.         exit(1);
  121.       }
  122.     } else if (!access("vox.bin.Z", R_OK)) {
  123.       if ((file = popen("zcat vox.bin.Z", "r")) == NULL) {
  124.         fprintf(stderr, "cannot popen input file vox.bin.Z (missing zcat?)\n");
  125.         exit(1);
  126.       }
  127.     } else {
  128.       fprintf(stderr, "cannot find vox.bin, vox.bin.gz, or vox.bin.Z\n");
  129.       exit(1);
  130.     }
  131.     using_pipe = 1;
  132. #else
  133.     fprintf(stderr, "cannot find vox.bin\n");
  134.     exit(1);
  135. #endif
  136.   } else {
  137.     using_pipe = 0;
  138.   }
  139.   vol_width = 128;      /* hard coded for demo */
  140.   vol_height = 128;
  141.   vol_depth = 64;
  142.   n_slices = 128;
  143.  
  144.   if (tex3dSupported) {
  145.     unsigned long size = (unsigned long) (vol_width * vol_height * vol_depth);
  146.  
  147.     vptr = (unsigned char *) malloc(size);
  148.  
  149.     fread(vptr, sizeof(char), size, file);
  150.     if (using_pipe) {
  151.       pclose(file);
  152.     } else {
  153.       fclose(file);
  154.     }
  155.  
  156.     /* size of voxels is twice as the size of vptr, to duplicate alpha value
  157.        = intensity */
  158.     voxels = (unsigned char *) malloc(2 * size);
  159.  
  160.     /* for now duplicate, alpha value = intensity */
  161.     vx = voxels;
  162.     vp = vptr;
  163.     for (i = 0; i < vol_width; i++)
  164.       for (j = 0; j < vol_height; j++)
  165.         for (k = 0; k < vol_depth; k++) {
  166.           *vx++ = *vp;
  167.           *vx++ = *vp++;
  168.         }
  169.  
  170.     free(vptr);
  171.   }
  172.   /* compute bounding box extents */
  173.  
  174.   bminx = -(float) vol_width / 2.0;
  175.   bmaxx = (float) vol_width / 2.0;
  176.   bminy = -(float) vol_height / 2.0;
  177.   bmaxy = (float) vol_height / 2.0;
  178.   bminz = -(float) vol_depth / 2.0;
  179.   bmaxz = (float) vol_depth / 2.0;
  180.  
  181.   bdiag = sqrt((bmaxx) * (bmaxx) + (bmaxy) * (bmaxy) + (bmaxz) * (bmaxz));
  182.  
  183.   /* compute view volume extents */
  184.  
  185.   left = -1.1 * bdiag;
  186.   right = 1.1 * bdiag;
  187.   bottom = -1.1 * bdiag;
  188.   top = 1.1 * bdiag;
  189.   nnear = -1.0 * bdiag;
  190.   ffar = 2 * 1.1 * bdiag;
  191.  
  192.   /* define the polygon dimensions on which the texture will be mapped */
  193.   /* xy plane */
  194.   slice_poly[0][0][0] = -vol_width / 2.0;
  195.   slice_poly[0][0][1] = -vol_height / 2.0;
  196.   slice_poly[0][0][2] = 0.0;
  197.  
  198.   slice_poly[0][1][0] = vol_width / 2.0;
  199.   slice_poly[0][1][1] = -vol_height / 2.0;
  200.   slice_poly[0][1][2] = 0.0;
  201.  
  202.   slice_poly[0][2][0] = vol_width / 2.0;
  203.   slice_poly[0][2][1] = vol_height / 2.0;
  204.   slice_poly[0][2][2] = 0.0;
  205.  
  206.   slice_poly[0][3][0] = -vol_width / 2.0;
  207.   slice_poly[0][3][1] = vol_height / 2.0;
  208.   slice_poly[0][3][2] = 0.0;
  209.  
  210.   /* yz plane */
  211.   slice_poly[1][0][0] = 0.0;
  212.   slice_poly[1][0][1] = -vol_height / 2.0;
  213.   slice_poly[1][0][2] = -vol_depth / 2.0;
  214.  
  215.   slice_poly[1][1][0] = 0.0;
  216.   slice_poly[1][1][1] = vol_height / 2.0;
  217.   slice_poly[1][1][2] = -vol_depth / 2.0;
  218.  
  219.   slice_poly[1][2][0] = 0.0;
  220.   slice_poly[1][2][1] = vol_height / 2.0;
  221.   slice_poly[1][2][2] = vol_depth / 2.0;
  222.  
  223.   slice_poly[1][3][0] = 0.0;
  224.   slice_poly[1][3][1] = -vol_height / 2.0;
  225.   slice_poly[1][3][2] = vol_depth / 2.0;
  226.  
  227.   /* zx plane */
  228.   slice_poly[2][0][0] = -vol_width / 2.0;
  229.   slice_poly[2][0][1] = 0.0;
  230.   slice_poly[2][0][2] = -vol_depth / 2.0;
  231.  
  232.   slice_poly[2][1][0] = -vol_width / 2.0;
  233.   slice_poly[2][1][1] = 0.0;
  234.   slice_poly[2][1][2] = vol_depth / 2.0;
  235.  
  236.   slice_poly[2][2][0] = vol_width / 2.0;
  237.   slice_poly[2][2][1] = 0.0;
  238.   slice_poly[2][2][2] = vol_depth / 2.0;
  239.  
  240.   slice_poly[2][3][0] = vol_width / 2.0;
  241.   slice_poly[2][3][1] = 0.0;
  242.   slice_poly[2][3][2] = -vol_depth / 2.0;
  243.  
  244.   /* texture coordinates */
  245.  
  246.   slice_tcoord[0][0][0] = 0.0;
  247.   slice_tcoord[0][0][1] = 0.0;
  248.   slice_tcoord[0][1][0] = 1.0;
  249.   slice_tcoord[0][1][1] = 0.0;
  250.   slice_tcoord[0][2][0] = 1.0;
  251.   slice_tcoord[0][2][1] = 1.0;
  252.   slice_tcoord[0][3][0] = 0.0;
  253.   slice_tcoord[0][3][1] = 1.0;
  254.  
  255.   slice_tcoord[1][0][1] = 0.0;
  256.   slice_tcoord[1][0][2] = 1.0;
  257.   slice_tcoord[1][1][1] = 1.0;
  258.   slice_tcoord[1][1][2] = 1.0;
  259.   slice_tcoord[1][2][1] = 1.0;
  260.   slice_tcoord[1][2][2] = 0.0;
  261.   slice_tcoord[1][3][1] = 0.0;
  262.   slice_tcoord[1][3][2] = 0.0;
  263.  
  264.   slice_tcoord[2][0][2] = 0.0;
  265.   slice_tcoord[2][0][0] = 0.0;
  266.   slice_tcoord[2][1][2] = 1.0;
  267.   slice_tcoord[2][1][0] = 0.0;
  268.   slice_tcoord[2][2][2] = 1.0;
  269.   slice_tcoord[2][2][0] = 1.0;
  270.   slice_tcoord[2][3][2] = 0.0;
  271.   slice_tcoord[2][3][0] = 1.0;
  272.  
  273. }
  274.  
  275. void
  276. randomTick(void)
  277. {
  278.   static unsigned int seed = 0;
  279.   static int changeSeed = 25;
  280.   float fltran;
  281.  
  282.   if (changeSeed++ >= 25) {
  283.     seed++;
  284.     if (seed > 256)
  285.       seed = 0;
  286.     changeSeed = 0;
  287.   }
  288.   fltran = (float) (rand_r(&seed) / 30000.0);
  289.  
  290.   anglex = (anglex > 360.0) ? 0.0 : (anglex + fltran);
  291.   angley = (angley > 360.0) ? 0.0 : (angley + fltran);
  292.   anglez = (anglez > 360.0) ? 0.0 : (anglez + fltran);
  293. }
  294.  
  295. void
  296. animate(void)
  297. {
  298.   randomTick();
  299.   glutPostRedisplay();
  300. }
  301.  
  302. void
  303. printCheatSheet(void)
  304. {
  305.   printf("\n\n-------------------------\n");
  306.   printf("OpenGL 3d texture example\n\n");
  307.  
  308.   printf("Keyboard shortcuts\n");
  309.   printf("s key   : zoom out (small)\n");
  310.   printf("l key   : zoom in  (large)\n");
  311.   printf("x key   : rotate about screen x\n");
  312.   printf("y key   : rotate about screen y\n");
  313.   printf("z key   : rotate about screen z\n");
  314.   printf("esc key : quit\n");
  315. }
  316.  
  317. void
  318. menu(int choice)
  319. {
  320.   /* simple GLUT popup menu stuff */
  321.   switch (choice) {
  322.   case SPIN_ON:
  323.     glutChangeToMenuEntry(1, "Random Spin OFF", SPIN_OFF);
  324.     glutIdleFunc(animate);
  325.     break;
  326.   case SPIN_OFF:
  327.     glutChangeToMenuEntry(1, "Random Spin ON", SPIN_ON);
  328.     glutIdleFunc(NULL);
  329.     break;
  330.   case MENU_HELP:
  331.     printCheatSheet();
  332.     break;
  333.   case MENU_EXIT:
  334.     exit(0);
  335.     break;
  336.   }
  337. }
  338.  
  339. void
  340. setMatrix(void)
  341. {
  342.   /* feel like using ortho projection */
  343.   glMatrixMode(GL_PROJECTION);
  344.   glLoadIdentity();
  345.   glOrtho(left, right, bottom, top, nnear, ffar);
  346.  
  347.   /* boring view matrix */
  348.   glMatrixMode(GL_MODELVIEW);
  349.   glLoadIdentity();
  350. }
  351.  
  352. void
  353. load3DTexture(void)
  354. {
  355.   printf("setting up 3d textures...\n");
  356.  
  357.   glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
  358.  
  359.   glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
  360.  
  361. #ifdef GL_EXT_texture3D
  362.   glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  363.   glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  364.   glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_WRAP_S, GL_REPEAT);
  365.   glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_WRAP_T, GL_REPEAT);
  366.   glTexParameteri(GL_TEXTURE_3D_EXT, GL_TEXTURE_WRAP_R_EXT, GL_REPEAT);
  367.  
  368.   glTexImage3DEXT(GL_TEXTURE_3D_EXT, 0, GL_LUMINANCE8_ALPHA8_EXT,
  369.     vol_width, vol_height, vol_depth,
  370.     0, GL_LUMINANCE_ALPHA, GL_UNSIGNED_BYTE, voxels);
  371. #endif
  372.  
  373.   if (tex3dSupported) {
  374.     /* enable texturing, blending etc */
  375.     glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  376.     glEnable(GL_BLEND);
  377.   }
  378.   setMatrix();
  379. }
  380.  
  381. void
  382. init(void)
  383. {
  384.   /* angle of rotation about coordinate axes */
  385.   anglex = angley = anglez = 0.0;
  386. }
  387.  
  388. void
  389. invert4d(float from[4][4], float to[4][4])
  390. {
  391.  
  392.   /* 4x4 matrix inversion routine */
  393.  
  394.   float wtemp[4][8];
  395.   register float m0, m1, m2, m3, s;
  396.   register float *r0, *r1, *r2, *r3, *rtemp;
  397.  
  398.   r0 = wtemp[0];
  399.   r1 = wtemp[1];
  400.   r2 = wtemp[2];
  401.   r3 = wtemp[3];
  402.   r0[0] = from[0][0];   /* build up [A][I] */
  403.   r0[1] = from[0][1];
  404.   r0[2] = from[0][2];
  405.   r0[3] = from[0][3];
  406.   r0[4] = 1.0;
  407.   r0[5] = 0.0;
  408.   r0[6] = 0.0;
  409.   r0[7] = 0.0;
  410.   r1[0] = from[1][0];
  411.   r1[1] = from[1][1];
  412.   r1[2] = from[1][2];
  413.   r1[3] = from[1][3];
  414.   r1[4] = 0.0;
  415.   r1[5] = 1.0;
  416.   r1[6] = 0.0;
  417.   r1[7] = 0.0;
  418.   r2[0] = from[2][0];
  419.   r2[1] = from[2][1];
  420.   r2[2] = from[2][2];
  421.   r2[3] = from[2][3];
  422.   r2[4] = 0.0;
  423.   r2[5] = 0.0;
  424.   r2[6] = 1.0;
  425.   r2[7] = 0.0;
  426.   r3[0] = from[3][0];
  427.   r3[1] = from[3][1];
  428.   r3[2] = from[3][2];
  429.   r3[3] = from[3][3];
  430.   r3[4] = 0.0;
  431.   r3[5] = 0.0;
  432.   r3[6] = 0.0;
  433.   r3[7] = 1.0;
  434.   if (r0[0] == 0.0) {   /* swap rows if needed */
  435.     if (r1[0] == 0.0) {
  436.       if (r2[0] == 0.0) {
  437.         if (r3[0] == 0.0)
  438.           goto singular;
  439.         rtemp = r0;
  440.         r0 = r3;
  441.         r3 = rtemp;
  442.       } else {
  443.         rtemp = r0;
  444.         r0 = r2;
  445.         r2 = rtemp;
  446.       }
  447.     } else {
  448.       rtemp = r0;
  449.       r0 = r1;
  450.       r1 = rtemp;
  451.     }
  452.   }
  453.   m1 = r1[0] / r0[0];   /* eliminate first variable */
  454.   m2 = r2[0] / r0[0];
  455.   m3 = r3[0] / r0[0];
  456.   s = r0[1];
  457.   r1[1] = r1[1] - m1 * s;
  458.   r2[1] = r2[1] - m2 * s;
  459.   r3[1] = r3[1] - m3 * s;
  460.   s = r0[2];
  461.   r1[2] = r1[2] - m1 * s;
  462.   r2[2] = r2[2] - m2 * s;
  463.   r3[2] = r3[2] - m3 * s;
  464.   s = r0[3];
  465.   r1[3] = r1[3] - m1 * s;
  466.   r2[3] = r2[3] - m2 * s;
  467.   r3[3] = r3[3] - m3 * s;
  468.   s = r0[4];
  469.   if (s != 0.0) {
  470.     r1[4] = r1[4] - m1 * s;
  471.     r2[4] = r2[4] - m2 * s;
  472.     r3[4] = r3[4] - m3 * s;
  473.   }
  474.   s = r0[5];
  475.   if (s != 0.0) {
  476.     r1[5] = r1[5] - m1 * s;
  477.     r2[5] = r2[5] - m2 * s;
  478.     r3[5] = r3[5] - m3 * s;
  479.   }
  480.   s = r0[6];
  481.   if (s != 0.0) {
  482.     r1[6] = r1[6] - m1 * s;
  483.     r2[6] = r2[6] - m2 * s;
  484.     r3[6] = r3[6] - m3 * s;
  485.   }
  486.   s = r0[7];
  487.   if (s != 0.0) {
  488.     r1[7] = r1[7] - m1 * s;
  489.     r2[7] = r2[7] - m2 * s;
  490.     r3[7] = r3[7] - m3 * s;
  491.   }
  492.   if (r1[1] == 0.0) {   /* swap rows if needed */
  493.     if (r2[1] == 0.0) {
  494.       if (r3[1] == 0.0)
  495.         goto singular;
  496.       rtemp = r1;
  497.       r1 = r3;
  498.       r3 = rtemp;
  499.     } else {
  500.       rtemp = r1;
  501.       r1 = r2;
  502.       r2 = rtemp;
  503.     }
  504.   }
  505.   m2 = r2[1] / r1[1];   /* eliminate second variable */
  506.   m3 = r3[1] / r1[1];
  507.   r2[2] = r2[2] - m2 * r1[2];
  508.   r3[2] = r3[2] - m3 * r1[2];
  509.   r3[3] = r3[3] - m3 * r1[3];
  510.   r2[3] = r2[3] - m2 * r1[3];
  511.   s = r1[4];
  512.   if (s != 0.0) {
  513.     r2[4] = r2[4] - m2 * s;
  514.     r3[4] = r3[4] - m3 * s;
  515.   }
  516.   s = r1[5];
  517.   if (s != 0.0) {
  518.     r2[5] = r2[5] - m2 * s;
  519.     r3[5] = r3[5] - m3 * s;
  520.   }
  521.   s = r1[6];
  522.   if (s != 0.0) {
  523.     r2[6] = r2[6] - m2 * s;
  524.     r3[6] = r3[6] - m3 * s;
  525.   }
  526.   s = r1[7];
  527.   if (s != 0.0) {
  528.     r2[7] = r2[7] - m2 * s;
  529.     r3[7] = r3[7] - m3 * s;
  530.   }
  531.   if (r2[2] == 0.0) {   /* swap last 2 rows if needed */
  532.     if (r3[2] == 0.0)
  533.       goto singular;
  534.     rtemp = r2;
  535.     r2 = r3;
  536.     r3 = rtemp;
  537.   }
  538.   m3 = r3[2] / r2[2];   /* eliminate third variable */
  539.   r3[3] = r3[3] - m3 * r2[3];
  540.   r3[4] = r3[4] - m3 * r2[4];
  541.   r3[5] = r3[5] - m3 * r2[5];
  542.   r3[6] = r3[6] - m3 * r2[6];
  543.   r3[7] = r3[7] - m3 * r2[7];
  544.   if (r3[3] == 0.0)
  545.     goto singular;
  546.   s = 1.0 / r3[3];      /* now back substitute row 3 */
  547.   r3[4] = r3[4] * s;
  548.   r3[5] = r3[5] * s;
  549.   r3[6] = r3[6] * s;
  550.   r3[7] = r3[7] * s;
  551.   m2 = r2[3];           /* now back substitute row 2 */
  552.   s = 1.0 / r2[2];
  553.   r2[4] = s * (r2[4] - r3[4] * m2);
  554.   r2[5] = s * (r2[5] - r3[5] * m2);
  555.   r2[6] = s * (r2[6] - r3[6] * m2);
  556.   r2[7] = s * (r2[7] - r3[7] * m2);
  557.   m1 = r1[3];
  558.   r1[4] = (r1[4] - r3[4] * m1);
  559.   r1[5] = (r1[5] - r3[5] * m1);
  560.   r1[6] = (r1[6] - r3[6] * m1);
  561.   r1[7] = (r1[7] - r3[7] * m1);
  562.   m0 = r0[3];
  563.   r0[4] = (r0[4] - r3[4] * m0);
  564.   r0[5] = (r0[5] - r3[5] * m0);
  565.   r0[6] = (r0[6] - r3[6] * m0);
  566.   r0[7] = (r0[7] - r3[7] * m0);
  567.   m1 = r1[2];           /* now back substitute row 1 */
  568.   s = 1.0 / r1[1];
  569.   r1[4] = s * (r1[4] - r2[4] * m1);
  570.   r1[5] = s * (r1[5] - r2[5] * m1);
  571.   r1[6] = s * (r1[6] - r2[6] * m1);
  572.   r1[7] = s * (r1[7] - r2[7] * m1);
  573.   m0 = r0[2];
  574.   r0[4] = (r0[4] - r2[4] * m0);
  575.   r0[5] = (r0[5] - r2[5] * m0);
  576.   r0[6] = (r0[6] - r2[6] * m0);
  577.   r0[7] = (r0[7] - r2[7] * m0);
  578.   m0 = r0[1];           /* now back substitute row 0 */
  579.   s = 1.0 / r0[0];
  580.   r0[4] = s * (r0[4] - r1[4] * m0);
  581.   r0[5] = s * (r0[5] - r1[5] * m0);
  582.   r0[6] = s * (r0[6] - r1[6] * m0);
  583.   r0[7] = s * (r0[7] - r1[7] * m0);
  584.   to[0][0] = r0[4];     /* copy results back */
  585.   to[0][1] = r0[5];
  586.   to[0][2] = r0[6];
  587.   to[0][3] = r0[7];
  588.   to[1][0] = r1[4];
  589.   to[1][1] = r1[5];
  590.   to[1][2] = r1[6];
  591.   to[1][3] = r1[7];
  592.   to[2][0] = r2[4];
  593.   to[2][1] = r2[5];
  594.   to[2][2] = r2[6];
  595.   to[2][3] = r2[7];
  596.   to[3][0] = r3[4];
  597.   to[3][1] = r3[5];
  598.   to[3][2] = r3[6];
  599.   to[3][3] = r3[7];
  600.   return;
  601. singular:
  602.   printf("ERROR : non_invertable transform\n");
  603.   return;
  604. }
  605.  
  606. void
  607. normalize(float *xn, float *yn, float *zn)
  608. {
  609.   double denom;
  610.  
  611.   denom = sqrt((double) ((*xn * *xn) + (*yn * *yn) + (*zn * *zn)));
  612.  
  613.   *xn = *xn / denom;
  614.   *yn = *yn / denom;
  615.   *zn = *zn / denom;
  616. }
  617.  
  618. int
  619. getViewAxis(void)
  620. {
  621.   float viewDir[3];
  622.   float mat[4][4], vinv[4][4];
  623.   float maxf, xy, yz, zx;
  624.   int im;
  625.  
  626.   /* out of 3 orthogonal set of planes in world coords,  find out which one
  627.      has maximum angle from the line of  sight.
  628.  
  629.      we will use these set of planes for creating polygonal slices thru the
  630.      voxel space */
  631.  
  632.   glGetFloatv(GL_MODELVIEW_MATRIX, (GLfloat *) mat);
  633.   invert4d(mat, vinv);
  634.  
  635.   viewDir[0] = -vinv[2][0];
  636.   viewDir[1] = -vinv[2][1];
  637.   viewDir[2] = -vinv[2][2];
  638.  
  639.   normalize(&viewDir[0], &viewDir[1], &viewDir[2]);
  640.  
  641.   xy = viewDir[2];      /* simplified because 0*xx + 0*yy + 1*zz */
  642.   yz = viewDir[0];
  643.   zx = viewDir[1];
  644.  
  645.   maxf = ABS(xy);
  646.   im = (xy < 0.0) ? _XY : _MXY;
  647.  
  648.   if (maxf <= ABS(yz)) {
  649.     maxf = ABS(yz);
  650.     im = (yz < 0.0) ? _YZ : _MYZ;
  651.   }
  652.   if (maxf <= ABS(zx)) {
  653.     maxf = ABS(zx);
  654.     im = (zx < 0.0) ? _ZX : _MZX;
  655.   }
  656.   return (im);
  657. }
  658.  
  659. #define DRAW_SLICE \
  660.             if (tex3dSupported) {\
  661.               glBegin(GL_POLYGON);\
  662.           for (p=0;p<4;p++) {\
  663.         glTexCoord3fv(slice_tcoord[myaxis][p]);\
  664.         glVertex3fv(slice_poly[myaxis][p]);\
  665.           }\
  666.           glEnd();\
  667.         } \
  668.             else {\
  669.               glBegin(GL_LINE_LOOP);\
  670.           for (p=0;p<4;p++){\
  671.         glVertex3fv(slice_poly[myaxis][p]);\
  672.           }\
  673.           glEnd();\
  674.         }
  675.  
  676. void
  677. drawScene(void)
  678. {
  679.  
  680.   int i, p;
  681.   float tc;
  682.   int viewAxis, myaxis;
  683.   int sign;
  684.  
  685.   static int myaxis_lut[] =
  686.   {0, 0, 1, 2, 0, 1, 2};
  687.  
  688.   /* clear background, z buffer etc */
  689.   glClear(GL_COLOR_BUFFER_BIT);
  690.  
  691.   glPushMatrix();
  692.  
  693.   /* apply all the modeling transformations */
  694.   glTranslatef(0.0, 0.0, -bdiag);
  695.   glRotatef(anglex, 1.0, 0.0, 0.0);
  696.   glRotatef(angley, 0.0, 1.0, 0.0);
  697.   glRotatef(anglez, 0.0, 0.0, 1.0);
  698.  
  699.   /* getViewAxis(), determines which set of polygons need to be used for
  700.      texturing, depending upon the viewing direction */
  701.  
  702.   viewAxis = getViewAxis();
  703.   myaxis = myaxis_lut[viewAxis];
  704.  
  705.   glColor3f(1.0, 1.0, 1.0);
  706.  
  707. #ifdef GL_EXT_texture3D
  708.   if (tex3dSupported)
  709.     glEnable(GL_TEXTURE_3D_EXT);
  710. #endif
  711.  
  712.   switch (viewAxis) {
  713.   case _XY:
  714.   case _MXY:
  715.  
  716.     sign = (viewAxis == _XY) ? 1 : -1;
  717.     tc = (viewAxis == _XY) ? 0.0 : 1.0;
  718.  
  719.     glTranslatef(0.0, 0.0, -sign * vol_depth / 2.0);
  720.     for (i = 0; i < n_slices; i++) {
  721.  
  722.       slice_tcoord[0][0][2] = slice_tcoord[0][1][2] =
  723.         slice_tcoord[0][2][2] = slice_tcoord[0][3][2] = tc;
  724.  
  725.       tc += sign * 1.0 / n_slices;
  726.       glTranslatef(0.0, 0.0, sign * vol_depth / (n_slices + 1.0));
  727.  
  728.       DRAW_SLICE;
  729.     }
  730.     break;
  731.   case _YZ:
  732.   case _MYZ:
  733.  
  734.     sign = (viewAxis == _YZ) ? 1 : -1;
  735.     tc = (viewAxis == _YZ) ? 0.0 : 1.0;
  736.  
  737.     glTranslatef(-sign * vol_width / 2.0, 0.0, 0.0);
  738.     for (i = 0; i < n_slices; i++) {
  739.  
  740.       slice_tcoord[1][0][0] = slice_tcoord[1][1][0] =
  741.         slice_tcoord[1][2][0] = slice_tcoord[1][3][0] = tc;
  742.  
  743.       tc += sign * 1.0 / n_slices;
  744.       glTranslatef(sign * vol_width / (n_slices + 1.0), 0.0, 0.0);
  745.  
  746.       DRAW_SLICE;
  747.     }
  748.     break;
  749.   case _ZX:
  750.   case _MZX:
  751.  
  752.     sign = (viewAxis == _ZX) ? 1 : -1;
  753.     tc = (viewAxis == _ZX) ? 0.0 : 1.0;
  754.  
  755.     glTranslatef(0.0, -sign * vol_height / 2.0, 0.0);
  756.     for (i = 0; i < n_slices; i++) {
  757.  
  758.       slice_tcoord[2][0][1] = slice_tcoord[2][1][1] =
  759.         slice_tcoord[2][2][1] = slice_tcoord[2][3][1] = tc;
  760.  
  761.       tc += sign * 1.0 / n_slices;
  762.       glTranslatef(0.0, sign * vol_height / (n_slices + 1.0), 0.0);
  763.  
  764.       DRAW_SLICE;
  765.  
  766.     }
  767.     break;
  768.   }
  769.  
  770. #ifdef GL_EXT_texture3D
  771.   if (tex3dSupported)
  772.     glDisable(GL_TEXTURE_3D_EXT);
  773. #endif
  774.  
  775.   glPopMatrix();
  776.  
  777.   glutSwapBuffers();
  778. }
  779.  
  780. void
  781. resize(int w, int h)
  782. {
  783.   /* things you do, when the user resizes the window */
  784.  
  785.   width = w;
  786.   height = h;
  787.  
  788.   glViewport(0, 0, w, h);
  789.   setMatrix();
  790. }
  791.  
  792. /* ARGSUSED1 */
  793. void
  794. keyboard(unsigned char c, int x, int y)
  795. {
  796.   /* handle key board input */
  797.  
  798.   switch (c) {
  799.   case 27:
  800.     exit(0);
  801.     break;
  802.   case 'x':
  803.     anglex += 1.0;
  804.     drawScene();
  805.     break;
  806.   case 'y':
  807.     angley += 1.0;
  808.     drawScene();
  809.     break;
  810.   case 'z':
  811.     anglez += 1.0;
  812.     drawScene();
  813.     break;
  814.   case 'm':
  815.     n_slices++;
  816.     drawScene();
  817.     break;
  818.   case 'e':
  819.     n_slices--;
  820.     if (n_slices < 4)
  821.       n_slices = 4;
  822.     drawScene();
  823.     break;
  824.   case 's':
  825.     if (left < right + 10.0) {
  826.       left -= 1.0;
  827.       right += 1.0;
  828.       bottom -= 1.0;
  829.       top += 1.0;
  830.       setMatrix();
  831.       drawScene();
  832.     }
  833.     break;
  834.   case 'l':
  835.     left += 1.0;
  836.     right -= 1.0;
  837.     bottom += 1.0;
  838.     top -= 1.0;
  839.     setMatrix();
  840.     drawScene();
  841.     break;
  842.   default:
  843.     break;
  844.   }
  845. }
  846.  
  847. void
  848. main(int argc, char **argv)
  849. {
  850.   int i, mode = GLUT_DOUBLE;
  851.  
  852.   glutInit(&argc, argv);
  853.  
  854.   for (i = 1; i < argc; i++) {
  855.     if (!strcmp(argv[i], "-no3Dtex")) {
  856.       tex3dSupported = 0;
  857.     } else if (!strcmp(argv[i], "-sb")) {
  858.       mode = GLUT_SINGLE;
  859.     }
  860.   }
  861.  
  862.   /* let glut do all the X Stuff */
  863.   glutInitDisplayMode(GLUT_RGB | mode);
  864.   glutCreateWindow("Voxel Head");
  865.  
  866.   /* init our variables, etc */
  867.   init();
  868.  
  869.   /* read texture data from a file */
  870.   readVoxelData();
  871.  
  872.   /* set up OpenGL texturing */
  873.   if (tex3dSupported)
  874.     load3DTexture();
  875.  
  876.   /* register specific routines to glut */
  877.   glutDisplayFunc(drawScene);
  878.   glutReshapeFunc(resize);
  879.   glutKeyboardFunc(keyboard);
  880.  
  881.   /* create popup menu for glut */
  882.   glutCreateMenu(menu);
  883.  
  884.   glutAddMenuEntry("Random Spin ON", SPIN_ON);
  885.   glutAddMenuEntry("Help", MENU_HELP);
  886.   glutAddMenuEntry("Exit", MENU_EXIT);
  887.  
  888.   glutAttachMenu(GLUT_RIGHT_BUTTON);
  889.  
  890.   /* loop for ever */
  891.   glutMainLoop();
  892. }
  893.